Nous sommes deux étudiants habitués aux métros parisiens, surchargés tout autant que le trafic routier, comme dans d’autres grandes villes. Nous pouvons donc affimer que se déplacer en vélo présente beaucoup d’avantages…mais malheureusement, un inconvénient existe : parfois il n’y a plus de vélos disponibles dans une station désirée, ou alors plus de places pour les rendre.
Nous avons lu un article d’un enseignant chercheur en mécanique des fluides qui a créé justement un logiciel de prédiction, intégré dans l’application «La Bonne Station» de Keolis, destinée aux utilisateurs du VCub a Bordeaux (équivalent des velibs a Paris) et nous avons trouvé intéressant cet enjeux de prévision. Son objectif est de maximiser l’efficacité des opérations de régulation des vélos: Quelles stations approvisionner en priorité ? Quelles stations agrandir en priorité ? Où positionner de nouvelles stations ? Quelles stratégies adopter pour être le plus efficace possible ?
Pour notre projet d’analyse prédictive reliée au vélibs, nous allons
travailler sur le jeu de données London bike.
Nous l’avons trouvé sur le site Kaggle.
Il correspond a des données de locations de vélos dans la ville de
Londres.
| timestamp | cnt | t1 | t2 | hum | wind | weather | holiday | weekend | season |
|---|---|---|---|---|---|---|---|---|---|
| 2015-01-04 00:00:00 | 182 | 3.0 | 2.0 | 93.0 | 6.0 | 3 | 0 | 1 | 3 |
| 2015-01-04 01:00:00 | 138 | 3.0 | 2.5 | 93.0 | 5.0 | 1 | 0 | 1 | 3 |
| 2015-01-04 02:00:00 | 134 | 2.5 | 2.5 | 96.5 | 0.0 | 1 | 0 | 1 | 3 |
| 2015-01-04 03:00:00 | 72 | 2.0 | 2.0 | 100.0 | 0.0 | 1 | 0 | 1 | 3 |
| 2015-01-04 04:00:00 | 47 | 2.0 | 0.0 | 93.0 | 6.5 | 1 | 0 | 1 | 3 |
Il contient \(17414\) observations, recensées sur deux ans, pendant chaque heure entre le 04/01/2015 à 00h00 et le 03/01/2017 à 23h00, et \(10\) variables :
\(\underline{Problématique}\)
Nous commencerons par chercher a prédire le nombre d’emprunts de vélos
des 4 derniers mois du jeu de données, avec des modèles qui utilisent
juste les variables explicatives.
Puis nous ferons un modele prédictif, afin de prevoir les futurs
emprunts de vélos.
Notre variable d’intérêt sera donc le nombre de vélos empruntés
«cnt».
Et nos variables explicatives seront les autres.
\(\underline{Traitement~des~données}\)
Nous commençons par verifier que nous n’avons pas de données
manquantes:
C’est bien le cas, donc nous n’avons pas besoin de faire d’imputation de celles-ci!
Ensuite, nous ajoutons plusieurs variables qui nous seront utiles pour la suite:
| timestamp | cnt | t1 | t2 | hum | wind | weather | holiday | weekend | season | day_of_month | month | day_of_week | hour | Posy | date | Year |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2015-01-04 00:00:00 | 182 | 3.0 | 2.0 | 93.0 | 6.0 | 3 | 0 | 1 | 3 | 4 | 1 | 1 | 0 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 01:00:00 | 138 | 3.0 | 2.5 | 93.0 | 5.0 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 1 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 02:00:00 | 134 | 2.5 | 2.5 | 96.5 | 0.0 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 2 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 03:00:00 | 72 | 2.0 | 2.0 | 100.0 | 0.0 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 3 | 0.0082418 | 2015-01-04 | 2015 |
| 2015-01-04 04:00:00 | 47 | 2.0 | 0.0 | 93.0 | 6.5 | 1 | 0 | 1 | 3 | 4 | 1 | 1 | 4 | 0.0082418 | 2015-01-04 | 2015 |
Voici le nombre de vélos empruntés au cours du temps:
Nos données ne recouvrent que deux ans, mais nous pouvons voir
notamment une saisonnalité.
Nous n’avons pas de changement signification d’une année à l’autre dans
la locations des vélos.
On remarque a priori qu’ils sont plus utilisés pendant l’été que
l’hiver.
Nous voyons deux valeurs anormales / deux pics élévés durant l’été 2015.
Les deux pics correspondent à des jours de grève des transports.
Regardons tout d’abord le nombre de vélos loués en fonction des différentes variables.
\(\underline{Variables~calendaires}\)
\(\bullet~Graphiques~annuel\)
Sur le premier graphique, nous voyons le nombre de locations de
vélos à chaque mois tout au long de l’année. La courbe noire correspond
a la premiere année de notre jeu de données (2015/2016), et celle en
bleue a la deuxieme année (2016/2017).
Nous voyons un pic élévé au milieu des deux années (7eme mois: juillet)
ce qui correspond aux étés. Tandis que les niveaux les plus bas sont en
hiver (1er mois: janvier) et sont inférieur à ceux des autres
saisons.
Cela nous confirme donc que les vélos sont moins empruntés pendant
l’hiver.
Sur le deuxieme graphique, nous avons la moyenne des deux années, et elle nous indique les même conclusions.
\(\bullet~Graphiques~hebdomadaire\)
Le premier graphique représente le nombre de locations de vélos
par jour(1=Dimanche…7=Lundi) en 1 semaine.
Nous remarquons que les vélos sont plus loués pendant la semaine que
pendant le week-end.
Nous avons également représenté ces resultats sous la forme d’un diagramme à barres.
\(\bullet~Graphiques~quotidien\)
Le premier graphique montre le nombre de locations de vélos par
heure pendant 1 jour.
La courbe bleue correspond aux jours de vacance, et celle en noire aux
autres jour.
Le deuxieme graphique nous montre le nombre de locations de vélos par
heure, pour chaque jours de la semaine.
Nous voyons clairement (et heureusement!!) que les locations pendant les non vacances, correspondent aux locations durant les jours de semaines. Et que les locations pendant les vacances, sont similaires a celles des jours de week end.
Par rapport aux vacances : nous remarquons qu’ils sont plus loués
aux alentours de 14heures.
Alors que le reste de l’année : ils sont plus loués aux moment des
heures de pointe (8 et 18 heures)
Et de même par rapport aux jours : En semaine, le nombre total de
locations est faible de 00h00 à 05h00 et augmente jusqu’au premier pic à
08h00, puis les déplacements se stabilisent et restent à peu près
constants jusqu’à 15h00.
Le deuxième pic est atteint à 17h00 et continue ensuite à baisser
jusqu’à environ 500 à 23h00.
Cela suggère potentiellement que les cyclistes utilisent les vélos
publics pour se rendre au travail vers 08h00 puis en revenir vers 17h00
en semaine.
Le week-end, les locations commencent à augmenter lorsque le soleil se
lève, et diminuent le soir.
\(\underline{Variables~liées~a~la~météo}\)
Nous remarquons que les vélos sont le plus loués lorsqu’il fait plus
chaud, et qu’il y a moins d’humidité.
La vitesse du vent n’a pas vraiment l’air d’influer sur les locations
par contre.
\(\underline{Corrélations}\)
Sur la matrice de corrélation, nous remarquons que notre variable d’interet cnt est le plus corrélée avec les variables t1, t2, hum et hour.
Sur le deuxieme “graphique” la proximité des points est
déterminée à l’aide du clustering multidimensionnel.
Les variables temperature réelle t1 et temperature ressentie
t2 apparaissent plus proche et sont jointes par un chemin plus
foncé, elles sont donc les plus fortement corrélées. Nous utiliserons
donc par la suite uniquement l’une des deux.
Les variables day_of_month, day_of_week,
holiday et weekend forment comme un petit groupe,
elles “donnent les mêmes explications”. Les autres variables
explicatives ne sont pas tellement corrélées entre elles.
\(\underline{Autocorrélations}\)
L’autocorrélation (ACF) d’une série fait référence au fait que
dans une série temporelle, la mesure d’un phénomène à un instant t peut
être corrélée aux mesures précédentes (au temps t-1, t-2, t-3, etc) ou
aux mesures suivantes (à t+1, t+2, t+3, …). L’autocorrélation mesure
essentiellement la similitude entre les observations en fonction du
décalage horaire entre elles. Une série autocorrélée est ainsi corrélée
à elle-même, avec un décalage (lag) donné.
Sur notre graphique ACF, l’autocorrélation est la plus forte pour les
décalages 1 et 24. (Et 0 mais ça ne signifie rien car l’autocorrélation
de la série avec elle-même, sans décalage, est forcément de 1!)
L’autocorrélation (PACF) partielle permet de mesurer
l’autocorrélation d’un signal pour un décalage k “indépendamment” des
autocorrélations pour les décalages inférieurs.
Sur notre graphique ACF, on a une forte autocorrélation pour un décalage
de 1, et de 24… Or l’autocorrélation pour le décalage de 24 pourrait
très bien être une conséquence de celle existant pour le décalage de 1,
c’est ce pourquoi nous regardons ensuite le graphique PACF. La valeur au
lag 1 est bien plus forte que pour les autres, ce qui nous laisse penser
que l’autocorrélation observée au décalage 24 était un effet résiduel de
l’autocorrélation pour un décalage de 1.
Si la série était stationnaire, nous verrions pratiquement toutes les lignes dans les intervalles de confiance a 95% représentés en verts. Or, nous voyons que chaque pic est hors de ces lignes.
Notre objectif est de développer des modèles basés sur l’ensemble de
données d’entraînement et de faire des prédictions du nombre horaire de
location de vélos à l’aide de l’ensemble de données de test.
Nous allons donc estimer des modèles sur des données d’entrainement et
faire des prévisions sur des données de test.
Nous commençons par divisé notre jeu de données en une partie pour
entrainer nos modèles, et une partie pour les valider.
Nous gardons les 4 derniers mois pour nos validations.
Nous avons 14488 observations pour nos tests, et 2926 pour nos validations.
Nous allons appliqué plusieurs modèles afin de prédire le nombre
total de locations de vélos par heure, pour les quatres dernier
mois.
Nous les évaluerons grâce a differentes valeurs :
Remarque:
- Nous avons choisi d’utiliser la MAPE et non la MAE (erreur absolue
moyenne), pour deux raisons:
La première est qu’un pourcentage est plus facile à comprendre, et la
seconde est que les scores MAPE peuvent être comparés entre différents
modèles. Cependant, si vos valeurs réelles auraient été proches de 0,
nous aurions choisis la MAE!!
- La MSLE (erreur logarithmique quadratique moyenne) réduit l’importance
des erreurs pour de grandes valeurs réelles, comme le fait deja la MAPE.
Elle ne s’exprime pas dans une unité simple, ce qui la rend difficile à
interpréter, ce pourquoi nous ne la choisissons pas non plus!
\(\underline{GLM}\)
RMSE = 1005.917
MAPE = 1.536
utilisateur système écoulé
0 0 0
Ce modèle est très rapide, mais le résultat n’est pas terrible..
\(\underline{SARIMA}\)
Series: bikeCount
ARIMA(5,1,0)
Coefficients:
ar1 ar2 ar3 ar4 ar5
0.3228 -0.3496 -0.0672 -0.0203 -0.0606
s.e. 0.0076 0.0080 0.0084 0.0079 0.0076
sigma^2 = 416609: log likelihood = -137366.9
AIC=274745.8 AICc=274745.8 BIC=274792.4
RMSE = 1328.121
MAPE = 4.406
utilisateur système écoulé
0.01 0.00 0.01
Le modèle SARIMA est une extension du modèle ARIMA.
Nous otenons une RMSE d’environ 1328 vélos, ce qui est encore moins bon
que pour notre modele de regression. On constate que le modèle est très
bon sur les données d’entraînement et moins bon sur les données de test.
On peut donc soupçonner des cas de sur-apprentissage.
\(\underline{Arbres~de~décision}\)
RMSE = 534.295
MAPE = 0.788
utilisateur système écoulé
0.24 0.02 0.25
Sur le premier graphique, nous pouvons visualiser l’arbre
obtenu.
Trois de nos variables hour, t1 et weekend
ont été utilisées, et 14 nœuds terminaux ont été générés.
Sur notre ensemble d’entrainement, nous voyons que la variable d’heure
hour est la première variable choisie pour la décision. Le
nombre de vélos empruntés attendu, lorsque hour<7 est de 197, et 29%
des échantillons tombent dans cette feuille. Lorsque hour est
inférieure à 7heure, moins de vélos sont utilisés, et lorsque
hour est >7h et <20h, la température commence à entrer en
jeu.
Sur le deuxieme graphique, nous pouvons voir l’erreur de
validation croisée (relative) pour chaque sous-arbre, du plus petit au
plus grand.
Dans “les coulisses”, rpart() applique automatiquement une validation
croisée pour élaguer l’arbre. L’axe x en bas correspond a la complexité
de l’arbre.
L’axe x en haut correspond au nombre de nœuds terminaux, donc la taille
de l’arbre.
L’axe y a l’erreur de convergence.
Un bon choix de cp pour l’élagage est souvent la valeur la plus à gauche
pour laquelle la moyenne se trouve sous la ligne horizontale. Nous
voyons que l’endroit optimal pour tailler l’arbre correspond a 12
nœuds.
Notre RMSE pour l’ensemble de test est d’environ 534
vélos!!
la MAPE vaut environ 0.79, et le temps d’execution est de 0.118
secondes.
Remarque: Les arbres de décision sont des algorithmes où les hyperparamètres sont essentiels en ce qui concerne leur performance. Par exemple, un nombre de profondeurs maximales (maxdepth) trop élevées peut entraîner un surajustement ou une variance élevée, et si ce nombre est trop faible cela peut entraîner a l’inverse un manque d’ajustement ou un biais élevé. Il faut donc faire attention car l’objectif est d’avoir les meilleurs resultats, mais l’arbre pourrait être suradapté.
\(\underline{GAM}\)
RMSE = 603.294
MAPE = 1.056
utilisateur système écoulé
0.04 0.00 0.05
Plus rapide (seulement 0.014 secondes), mais la RMSE est d’environ 603 vélos, et la MAPE d’environ 1.06 : c’est un peu moins bien que l’arbre !
\(\underline{Randomn~forest}\)
RMSE = 557.427
MAPE = 0.835
utilisateur système écoulé
0.16 0.11 0.14
Notre RMSE pour l’ensemble de test est d’environ 549 vélos, ce qui est un tout petit peu plus que pour l’abre de décision, mais ils ont la même valeur de MAPE, environ 0.79, et nos forets aleatoires sont plus rapide (0.028 secondes)
\(\underline{XGBoost}\)
RMSE = 372.815
MAPE = 0.718
utilisateur système écoulé
0.04 0.00 0.02
La RMSE vaut environ 373 vélos et la MAPE environ 0.718!!! Nous obtenons les meilleurs resultats, l’extreme gradient boosting dépasse tout le monde! De plus le temps d’execution est seulement de 0.002 secondes.
\(\underline{Perceptron}(MLP)\)
Entrainons le perceptron multi-couches en faisant varient le nombre de neurones, le nombre de couches afin de trouver un modèle qui prédit mieux.
La fonction RELU \(f(x)=x^{+}=max(0,x)\) sera notre fonction d’activation. On utilisera le package Keras.
128 neurones à la première couche et 1 à la dernière.
RMSE = 909.46
MAPE = 1.345
utilisateur système écoulé
0.44 0.82 0.32
500 neurones à la première et deuxième couche et 1 à la dernière.
RMSE = 669.07
MAPE = 0.542
utilisateur système écoulé
0.67 1.03 0.53
600 neurones à la première,deuxième et troisième couche et un à la dernière.
RMSE = 635.375
MAPE = 0.364
utilisateur système écoulé
1.19 1.05 0.59
On constate que le meilleur modèle est le 3 en termes de RMSE et MAPE. Donc plus le modèle aura des couches et de neurones mieux sera les résultats.
Nous allons maintenant faire un modele prédictif pour essayer de prevoir le jour même, les 24 points (1 point par heure) pour le lendemain, avec nos données allant jusqu’à la veille.
\(\underline{Prophet}\)
Prophet se base sur un modèle additif où les tendances non linéaires
correspondent aux saisons annuelles, hebdomadaires et quotidiennes,
ainsi qu’aux effets des vacances.
Le modèle fonctionne mieux avec les séries temporelles qui ont de forts
effets saisonniers et les données historiques de plusieurs
saisons.
Il examine les données antérieures afin de pouvoir analyser le
futur.
Il utilise un modèle de série temporelle qui peut être décrit par 3
composantes principales du modèle: les tendances, les saisons et les
vacances, selon l’équation suivante :
\(y(t) = g(t) + s(t) + h(t) + \epsilon(t)\)
\(g(t)\) : changements non périodiques des valeurs de la série \(s(t)\) : changement périodique (saisonnier, hebdomadaire, annuel) \(h(t)\) : l’effet des jours fériés sur le calendrier, pouvant créer des irrégularités \(\epsilon(t)\) : changements inhabituels non pris en compte par le modèle
Prophet utilise des régresseurs de température et d’humidité. Tout
cela nous semble très interessant dans notre cas.
Nous allons donc faire des prédictions pour le prochain jour pour chaque
heure (ou pour 24 heures).
ds yhat
17414 2017-01-03 23:00:00 158.7483
17415 2017-01-04 00:00:00 766.0758
17416 2017-01-04 01:00:00 1293.9043
17417 2017-01-04 02:00:00 1499.2879
17418 2017-01-04 03:00:00 301.4787
17419 2017-01-04 04:00:00 -115.5142
17420 2017-01-04 05:00:00 798.4119
17421 2017-01-04 06:00:00 1346.4293
17422 2017-01-04 07:00:00 1761.7837
17423 2017-01-04 08:00:00 3551.3491
17424 2017-01-04 09:00:00 2239.1037
17425 2017-01-04 10:00:00 1764.9763
17426 2017-01-04 11:00:00 1911.5235
17427 2017-01-04 12:00:00 1076.4066
17428 2017-01-04 13:00:00 2219.9535
17429 2017-01-04 14:00:00 1773.2393
17430 2017-01-04 15:00:00 1525.1560
17431 2017-01-04 16:00:00 2210.3772
17432 2017-01-04 17:00:00 3270.4694
17433 2017-01-04 18:00:00 3459.5608
17434 2017-01-04 19:00:00 2913.8282
17435 2017-01-04 20:00:00 2161.0340
17436 2017-01-04 21:00:00 1736.3712
17437 2017-01-04 22:00:00 1255.1763
17438 2017-01-04 23:00:00 635.5908
Nous retrouvons la prédiction de nos 24 points pour le lendemain.
Sur ce graphique les points en noirs representent les données réelles,
et les lignes bleues sont les valeurs prédites. Nous le dernier pic bleu
correspond a la prediction future de la journée suivante.
Nous voyons les performances du modèle grâce au nuage de points, entre
la valeur prédite et la valeur réelle.
La droite ajustée (en rouge) a une pente positive, ce qui signifie que
plus les valeurs “actual” sont grandes, plus les valeur “pred” le sont
aussi.
Adjusted R-squared: 0,6261
[R-squared est un nombre compris entre 0 et 1, et plus il est grand,
meilleur sera le modèle résultant] Notre résultat n’est pas le plus
optimal.
Nous terminons cette partie modelisation par un test ARMA avec
Prophet.
Nous allons essayer de faire une rediction du futur sur un an avec la
moyenne mobile.
À présent agrégions les différents modèles obtenu jusqu’à présent, à savoir l’arbre de décision, le modèle additif généralisé, la forêt aléatoire, le gradient boosting et le perceptron multi-couches. On utilisera le package opera pour l’agrégation..
[1] 2926 5
Aggregation rule: EWA
Loss function: squareloss
Gradient trick: FALSE
Coefficients:
rpart gam rf xgboost mlp
6.35e-216 0 5.02e-253 1 0
Number of experts: 5
Number of observations: 2926
Dimension of the data: 1
rmse mape
EWA 373 0.716
Uniform 433 0.518
On peut toute de suite constater que l’algorithme n’a utilisé que les meilleurs modèles. Le gradient boosting domine tous les autres experts. Le perceptron multi-couches apparaît tout au début en éclipsant les autres experts.
Le RMSE et le MAPE du modèle agrégé sont de 373 et 0.715. Soit le meilleur mape que ceux de tous nos modèles précédents et le deuxième meilleur rmse derrière xbgboost de rmse 372.
Nous remarquons que pour certains modèles, les prévisions du nombre de vélos prennent parfois des valeurs négatives, ce qui n’est pas possible mais comme ce sont des prévisions, cela fait partie de “l’erreur”
| RMSE | MAPE | Temps d’execution | |
|---|---|---|---|
| GLM | 1005.91704439951 | 1.53614602226818 | 0.07 |
| SARIMA | 1328.12067903949 | 4.40586219999758 | 0.02 |
| Arbre de décision | 534.294666099615 | 0.787693158684004 | 0.118 |
| Modèle Additif Généralisé | 603.293713663435 | 1.05554785178295 | 0.007 |
| Foret Aléatoire | 557.427418914906 | 0.835415619899165 | 0.029 |
| eXtreme Gradient Boosting | 372.81500435415 | 0.717676653923912 | 0.002 |
| MLP | 635.374503996781 | 0.363700834869851 | 0.67 |
| EWA | 373 | 0.717 | NA |
| Uniform | 575 | 0.498 | NA |
Nous obtenons des MAPE < 10%, ce qui est très bien !
\(\underline{Validation~du~modèle~final}\)
Nous pouvons constater que c’est avec le modèle XGBoost que nous obtenu
les meilleurs résultats, et c’est celle qui a également nécessité le
moins de temps de calcul, car elle demande moins de paramètres à former
que les autres modèles.
Les algorithmes XGBoost et Random Forest souffrent d’une variance élevée, alors on pourrait ajouter plus de données et augmenter l’effet de la régularisation… Nous avons entrainé nos modèles sur des données d’un an et 8 mois, et validé sur un ensemble de test de 4 mois. Les vélibs sont moins loués pendant l’hiver que pendant le reste de l’année. Peut-être que l’ajout de plus de données pendant l’hiver aidera à améliorer les performances pendant cette période, ou faire un sur-échantillonnage.
En conclusion générale nous sommes ravis d’avoir pu tester ces modèles sur un jeu concret de donnée.